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Detecting gravitationally lensed supernovae is among the biggest challenges 
inastronomy. It involves a combination of two very rare phenomena: 
catching the transient signal of astellar explosion ina distant galaxy and 
observing it through a nearly perfectly aligned foreground galaxy that 
deflects light towards the observer. Here we describe how high-cadence 
optical observations with the Zwicky Transient Facility, with its unparalleled 
large field of view, led to the detection ofa multiply imaged type la supernova, 
SN Zwicky, also known as SN 2022qmx. Magnified nearly 25-fold, the system 
was found thanks to the standard candle nature of type la supernovae. 
High-spatial-resolution imaging with the Keck telescope resolved four 
images of the supernova with very small angular separation, corresponding 


to an Einstein radius of only 0; = 0.167” and almost identical arrival times. The 
small 6, and faintness of the lensing galaxy are very unusual, highlighting the 
importance of supernovae to fully characterize the properties of galaxy-scale 


gravitational lenses, including the impact of galaxy substructures. 


Our understanding of gravitational lensing due to the curvature 
of spacetime, and the analogy with the deflection of light in optics, 
dates back to the work of Einstein in 1936’. In this pioneering work 
he considered the case where both the lens and the magnified back- 
ground source were stars in our Galaxy. Einstein concluded that the 
deflection angles were too small to be resolved with astronomical 
instruments. It was Zwicky’ who one year later pointed out that, if the 
source was extragalactic, entire galaxies or clusters of galaxies could 
be considered as gravitational deflectors. Hence, the image separation 


between multiple images of the source could be large enough to 
be resolved by astronomical facilities, as the size of the image sepa- 
ration scales with the lens mass and distance as the Einstein radius, 


l 1 1 
a nf M \2( Ds \ 2 Dis 2 ; 
Og x 0.9 (sc) (=) (2 ) , where M, is the mass of the Sun, M, 
and D, are the lensing mass and lens angular size distance and D, and 
D,, are the distances from the observer to the source and between the 
lens and the source, respectively. 


A full list of affiliations appears at the end of the paper. 


e-mail: ariel@fysik.su.se 


Nature Astronomy 


Article 


https://doi.org/10.1038/s41550-023-01981-3 


Restframe wavelength (A) 


3,000 4,000 5,000 6,000 
-14 L L L i 
P60/SEDM, Aug. 21.2 
NOT/ALFOSC, Aug. 22.9 
fe 5 4 
Ta 
N 
z SN 2012cg (shifted to z = 0.354) 
oO 
p 
2 
© 
N 
€ 
© 
+ 
us 16 4 
D 
2 
Keck/LRIS, Aug. 23.3 
-17 T T T T Vy T 
4,000 5,000 6,000 7,000 8,000 9,000 
Observed wavelength (A) 
[O] Ca ıı H&K Nai D Ha [N u] [S 11] 
us 
T Host 
oO 
a 
© Ui Ibo ii ha aa 
E 
[o] 
Z Lens PA 
T T T T T T T T T 
3,700 3,800 3,900 4,000 5,800 5,900 6,000 6,600 6,700 


Restframe wavelength (A) 


Fig. 1| Spectroscopic identification of SN Zwicky as an SN Ia and redshift 
measurements of the SN host galaxy and the intervening lensing galaxy. 
The SN spectra obtained with the P60, Nordic Optical Telescope (NOT) and Keck 
telescopes (black lines) are best fitted by anormal SN Ia spectral template. The 
green line shows a comparison with the nearby type la SN 2012cg” at a similar 
restframe phase, redshifted to z = 0.354. The SN flux peaked on 2022 August 17. 
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The bottom panels show a zoomed-in view of a VLT/MUSE spectrum from 2022 
September 30, displaying narrow absorption and emission lines, from which 
precise redshifts of the lens (z= 0.2262) and host (z= 0.3544) galaxies were 
determined. [O 11], Ca 11 H and K, NaID, Ha, [N 11] and [S 11] lines can be seen at the 
restframe of the lens (blue lines) and host (red lines) galaxies. ALFOSC, Alhambra 
Faint Object Spectrograph and Camera. 


Strongly lensed supernovae in the era of 
wide-field time-domain surveys 

Besides the many observations of lensed galaxies and quasars, the feasi- 
bility of observing strong gravitational lensing of explosive transients in 
the distant universe has only been demonstrated inrecent years (refs. 3-5 
and references therein). PS1-10afx was the first highly magnified type 
la supernova (SN la) discovered. However, the lensing interpretation 
was made three years after the explosion®”’; by then the SN was too faint 
to resolve multiple images. Since then, the use of wide-field optical 
cameras in robotic telescopes at the Palomar Observatory has led to 
notable breakthroughs. In ref. 8 we reported the discovery of a multiply 
imaged SN la, iPTFl6geu (SN 2016geu), by the intermediate Palomar 
Transient Factory (iPTF), a time-domain survey using a 7.3%? camera 
on the P48 (1.2 m) telescope from 2013 to 2017. In 2018, anew camera 
was installed’ with a field of view of 47°. The project, known as the 
Zwicky Transient Facility (ZTF)'°", has been monitoring the northern 
sky with a 2-3 d cadence in at least two optical filters for the past four 
years”. The very large sky coverage makes ZTF well suited to search 
for rare phenomena, such as gravitational lensing of supernovae. On 
the other hand, the distance (redshift) probed by ZTF is limited by the 


relatively small mirror of the telescope, light pollution, non-optimal 
atmospheric conditions and only having three optical filters at the 
P48 telescope. Furthermore, ZTF typically obtains an image quality 
(angular resolution) of 2” full-width at half-maximum (FWHM), and 
the camera has relatively large 1” pixels. Hence, in most instances, it is 
practically impossible to spatially resolve multiple-image systems with 
ZTF. Instead, the search for lensed sources makes use of the standard 
candle nature of type la supernovae, that is, they have nearly identi- 
cal peak luminosities. These explosions are used as accurate distance 
estimators in cosmology, which led to the discovery of the accelerated 
expansion of the universe (ref. 13 and references therein). 

In addition to an imaging survey telescope, ZTF has access to a 
low-spectral-resolution integral-field spectrograph, the SED Machine 
(SEDM)", on the neighbouring 1.5 m telescope at Palomar (P60), used 
to spectroscopically classify about ten supernovae every night as part 
of the Bright Transient Survey (BTS), where transients brighter than 
19 mag are classified within timescales of a few days, aiming to obtain 
>95% spectroscopic completeness to 18.5 mag or brighter’. Besides 
providing the classification of the transients, the SEDM spectrum is 
used to measure the SN redshift. 
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The discovery of SN Zwicky 

Lensed system candidates are selected by ZTF for further spectroscopic 
screening when an SN lais found at a redshift above z = 0.2, where there 
is essentially negligible sensitivity for detection inthe BTS, unless the 
SNis greatly magnified by an intervening deflector. This was the case 
for SN Zwicky (also known as ZTF22aaylInhq and SN 2022qmx), located 
at right ascension 17 h 35 min 44.32 s and declination 4° 49’ 56.90” 
(J2000.0), where an SEDM spectrum from 2022 August 21 showed it 
to bean SN la at z= 0.35, as shown in Fig. 1. At that point we alerted the 
SN community to the discovery of a lensed SN la". 

Spectroscopic observations following the time evolution of the SN 
were carried out using multiple facilities: the 2.56 m Nordic Optical Tel- 
escope in the Canary Islands, the Keck observatory in Hawaii, the 11 m 
Hobby-Eberly Telescope at the McDonald Observatory in Texas and 
the European Southern Observatory (ESO)’s 8 m Very Large Telescope 
(VLT) at the Paranal Observatory in Chile. In particular, multiple narrow 
emission and absorption lines of the SN host galaxy were found with the 
Low Resolution Imaging Spectrometer (LRIS)/Keck and the Multi-Unit 
Spectroscopic Explorer (MUSE)/VLT, refining the source redshift to 
z = 0.3544, as shown in the bottom panels of Fig. 1. As the SN faded 
and the foreground galaxy spectral energy distribution became more 
prominent, the Ca 11 doublet AA3,933, 3,968 was found in absorption 
lines redshifted to z = 0.22615, the location of the deflecting galaxy. 


Follow-up observations 
The discovery from ZTF wasalso followed up with high-spatial-resolution 
instruments. Observations with laser guide star enhanced seeing at the 
VLT with the High Acuity Wide-field K-band Imager (HAWK-I) imag- 
ing camera in the near infrared, and optical spectrophotometry with 
MUSE, reduced the point spread function (PSF) width to about 0.4”. 
However, this was still not enough to resolve the system. On 2022 Sep- 
tember 15, multiple images of the system were resolved at near-infrared 
wavelengths at the Keck telescope, using the Laser Guide Star Aided 
Adaptive Optics (LGSAO) with Near-IR Camera 2 (NIRC2)”, yielding an 
image quality of 0.09” FWHM intheJ band centred at 1.2 um, shown in 
Fig. 2, where the four SN images are labelled A-D. 

On September 21, following our announcement of the discovery", 
a previously approved programme aimed to target lensed supernovae 
by the LensWatch collaboration resolved the multiple images of SN 
Zwicky using the optical filters F475W, F625W and F814W (where the 
names correspond to the approximate location of the central wave- 
length in nanometres) on the UVIS/WFC3 camera on the Hubble Space 
Telescope (HST)"*. A detailed description of the HST observations of 
SN Zwicky is presented in ref. 19. 


Results 

Figure 3 shows the unresolved photometric ground-based observations 
collected at P48 and the Liverpool Telescope in the Canary Islands 
between 2022 August 1 and October 30. These were used to estimate the 
peak flux and lightcurve properties of the SN with the SALT2 lightcurve 
fitting tool”, including corrections for lightcurve shape and colour 
excess given the SN redshift, as well as the extinction in the Milky Way 
in the direction of the SN . Furthermore, the four resolved SN images 
were used to explore the possibility of additional reddening by dust 
inthe lensing galaxy. Unaccounted-for dimming of light would lead to 
an underestimation of the lensing amplification. The HST and NIRC2 
observations for each SN image were compared with the SN la spectral 
template from ref. 21, allowing for possible differential dust extinction 
in the lens following the reddening law in ref. 22. 

No evidence for differential extinction between the different 
images was found. The lightcurve fit model included the four indi- 
vidual SN images, described by the SALT2 model with arbitrary time 
delays, but otherwise sharing the same lightcurve shape and colour 
parameters, x, and c. The time delays were constrained by a prior onthe 
image flux ratios from the NIRC2 observations shown in Fig. 2. We find 
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Fig. 2 | Image of the field of SN Zwicky using pre-explosion g- and r-band 
images from ZTF. Top left: a2’ x 2’ section of the ZTF g- and r-band pre-SN 
images (FWHM 2.3”), centred on the location of SN Zwicky. Top right: a 
zoomed-in composite image of SN Zwicky using adaptive optics (AO)-enhanced 
VLT MUSE (g/r-band) and HAWK-I (J-band) observations on September 2 

and 4 (FWHM -0.4”). Bottom: a 2” x 2” portion of the Keck/NIRC2 LGSAO 

J-band image (FWHM 0.09”), resolving the four multiple images of SN Zwicky 
(labelled A, B, C, D). The blue dashed ellipse shows the critical line of the lens, 
corresponding to the inferred 6; = 0.167” (0.6 kpc at z = 0.2262), enclosing the 
lens M= (7.82 + 0.06) x 10° M,. The host galaxy nucleus is located 1.4” to the 
northeast of the lens, implying that SN Zwicky exploded at a projected distance of 
7 kpc from the centre of its host galaxy. 


Atay = —0.4 + 2.9, Atc = -0.1 + 2.3 and At,» = -0.1 + 2.7 (in units of days), 
where the indices A-D refer to the SN images in Fig. 2. The resulting 
lightcurve fit is shown in Fig. 3 and compared with the SALT2 model. 
The small time delays are also consistent with the spectra of SN Zwicky 
shown in Fig. 1 being at a single SN phase. 

The fitted SN model parameters are x, = 1.083 + 0.094 and 
c = -0.007 + 0.007. The lack of colour excess confirms that differen- 
tial extinction is negligible. Since the lightcurve parameter errors 
do not include the model covariance, we conservatively add the 
SALT2 model error floor of o(x,) = 0.1 and o(c) = 0.027 mag (ref. 23) in 
quadrature to the fit errors. Using the inferred apparent magnitude 
and the SALT2 parameters above, we find a total magnification of 
Am=—3.44 + 0.14 mag, assuming standard cosmological parameters 
from ref. 24 and a restframe B-band SN la peak absolute magnitude of 
-19.4 mag for the average SN la lightcurve width and colour, and intrin- 
sic brightness scatter of 0.1 mag. Since the inferred stellar mass of the 
host galaxy is M, < 10° M,, mass-step corrections for the SN Ia abso- 
lute magnitude as suggested in ref. 25 are not required. In summary, 
we find that, including the four images, SN Zwicky is 23.7 + 3.2 times 
brighter than the observed flux of normal type la supernovae at the 
same redshift, after applying colour and lightcurve shape corrections. 

The Keck NIRC2 J-band image was used to obtain a lens model 
to account for the observed SN image positions, irrespective of their 
fluxes. Assuming a singular isothermal ellipsoid”’” for the lens poten- 
tial, we report an ellipticity €, = 0.35 + 0.01 and 6, = 0.1670 + 0.0006”. 
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Fig. 3 | Multicolour lightcurve of SN Zwicky showing that the supernova is 
3.5 mag brighter than an unlensed SN at the same redshift. The magnitudes 
are measured with respect to time of maximum light (modified Julian Date 
59808.67, corresponding to 2022 August 17) in ZTF g and r and Liverpool 
Telescope (LT) g, r,i, z filters. The solid lines show the SALT2” model with the 
best fit to the spatially unresolved data. The blue dashed lines indicate the 
expected lightcurves at z= 0.354 (without lensing), where the bands represent 
thes.d. of the brightness distribution for type la supernovae. To fit the observed 
lightcurves, a brightness increase corresponding to 3.5 mag is required. Also 


shown, as dotted lines, are the modelled individual lightcurves for the four SN 
images A-D. The flux ratios were measured from the NIRC2 data in Fig. 2. From 
these lightcurves we extract the time delays between the images, all in units 

of days, At, = -0.4 + 2.9, At,c = -0.1 + 2.3 and At,» = -0.1 + 2.7, as described in 
Supplementary Information. The shaded areas in the lower panels indicate the 
regions with data outside the phase range where the SALT2 model is defined, 
and therefore excluded from the lightcurve fit analysis. The error bars 
correspond to1s.d. 


The mass enclosed within the ellipse with semi-major axis 0.78 kpc and 
semi-minor axis 0.51 kpc is M = (7.82 + 0.06) x 10° M,. The lens model 
predictions for the time delays are in excellent agreement with the fit- 
ted values from the lightcurves in Fig. 3. Further details regarding the 
lens modelling are presented in Methods. 

Interestingly, the individual image magnifications predicted for 
SN Zwicky by the smooth macro lens model are inconsistent with the 
observed flux ratios. According to the lensing model, the observed 
fluxes of the SN images A and Care factors of >4 and >2 too large, respec- 
tively, compared with images B and D. Given the small time delays, this 
discrepancy cannot be accounted by different phases between the 
SN images. Other explanations need to be considered: for example, 
excess magnification and demagnification from milli- and microlensing 
effects arising from stars and substructure in the lens galaxy**”’. Since 
microlensing effects are capable of perturbing magnifications without 
altering image locations, these were also put forward to explain the 
observed flux ratios of iPTFl6geu*’, displaying differences between 
the observed and modelled image flux ratios of similar magnitudes. 
Probing microlensing in these central regions opens a new window 
to directly measure the central stellar initial mass function” and test 
claims that the initial mass function may be heavier in the centres of 
galaxies”. As detailed in Methods, the lack of time dependence in the 
image flux ratios, or anomalous variation inthe unresolved lightcurves, 
gives a lower limit for the substructure masses of 0.02 M,, ifthe discrep- 
ancy fromthe smooth lens model is caused by microlensing. From the 


lack of further image splitting of the four individual SN images, we infer 
an upper limit for the substructure mass of 3 x 10’ M,. 

To check the impact of added macro lens model complexity, we 
have studied cases where the lens mass distribution is modelled with 
two matter components; one where the surface mass density follows 
the lens light distribution (with an arbitrary mass-to-light ratio, pos- 
sibly interpreted as a baryonic mass component), and a second one 
introducing a dark matter halo with additional flexibility on density 
profiles. In spite of the added extra complexity, the quality of the fit 
to the SN image positions does not improve, and induces shifts in the 
predicted flux ratios only below 5%, that is, very small in comparison 
with the observed flux ratio anomalies. 

The demonstrated ability to discover multiply imaged super- 
novae makes it feasible to accomplish Refsdal’s pioneering proposal 
3 to use time delays for strongly lensed supernovae to measure the 
Hubble constant. This will require systems with time delays of several 
days, that is, longer than for SN Zwicky. The small physical scale of the 
lens probed by SN Zwicky, as well as iPTFl6geu®, make these supernovae 
unique tracers for uncovering a population of systems that otherwise 
would remain undetected, as shown in Fig. 4, probing the mass distribu- 
tion of the central, densest, regions of lensing galaxies. Both of these 
multiply imaged type la supernovae were identified without prese- 
lections on, for example, association with bright galaxies or clusters, 
emphasizing the importance of untargeted surveys for unexpected 
discoveries. 
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Fig. 4 | Stellar mass versus Einstein radius for lens galaxies discovered 

in galaxy surveys, demonstrating that SN Zwicky, iPTF16geu and the 
unresolved lensed supernova PS1-10afx point to a poorly known population 
of small-image-separation lensing systems. Strongly lensed galaxy systems 
are represented by yellow diamonds for the BOSS Emission-Line Lens Survey 
(BELLS) sample”, blue triangles for Sloan Lens ACS (SLACS) lenses” and purple 
squares for SL2S lenses”. The green circles correspond to lensed quasars”, of 
which the filled circles have been detected optically and the open circles through 
radio emission. The shaded grey contours show the 90% and 68% confidence 
levels for the full sample of 155 lensed galaxies and 45 lensed quasars. The lensed 
supernova data are presented as median values + 1 s.d. For the unresolved lensed 
supernova PS1-10afx, only an upper limit of the Einstein radius is available’. The 
stellar mass derivation for SN Zwicky is detailed in Methods. 


Methods 

Supernova survey and follow-up 

The ZTF has been monitoring the transient sky at optical wave- 
lengths since 2018" ""**. SN Zwicky™ was discovered under the BTS 
programme’. The first detection of the SN was ina ZTF g-band image 
from 2022 August 1. It was saved to the BTS as an SN candidate by 
on-duty scanners on August 3” and subsequently assigned to the 
queue for spectroscopic follow-up with SEDM mounted on P60 under 
standard BTS protocols. The SEDM spectrum, obtained on August 
21°°, shows an excellent match to anormal SN Ia at a redshift of z= 0.35 
close to maximum light. The redshift and spectral classification were 
confirmed with a higher-resolution spectrum obtained at the Nordic 
Optical Telescope in La Palma on the following night. We followed up SN 
Zwicky with P48 inthe g andr bands. For our analysis we use the forced 
photometry provided by the Image Processing and Analysis Center as 
detailed in ref. 34. Observations in the Sloan Digital Sky Survey g, r,iand 
z filters were taken with the 10:0 optical imager on the Liverpool Tel- 
escope”. The Liverpool Telescope photometric data are processed with 
custom data-reduction and image-subtraction software (K. Taggart 
et al., manuscript in preparation). Image subtraction is performed 
using the Panoramic Survey Telescope and Rapid Response System 
1 (Pan-STARRSI) reference image. Images were stacked using SWarp 
to combine multiple exposures where required. The photometry is 
measured using a PSF fitting methodology relative to Pan-STARRS1 
standards and is based on techniques in ref. 38. 


Lightcurve fit, magnification and time-delay inference 

We used the publicly available, Python-based software SNTD” for 
inferring the restframe B-peak magnitude, the lightcurve shape and 
colour SN Ia SALT2 parameters” and the time delays between the SN 
images. Data points with >30 detections from the g and r filters from 
P48 andtheg, r, i,z filters from the Liverpool Telescope were included 
in the fit. We adopted an iterative procedure in two steps. First, all the 
data were used. Next, only data points in the range where the SALT2 


model is defined, that is, -20 to +50 d, were kept for the second itera- 
tion. The final lightcurve fit parameters were derived from data inthis 
phase range. 

We estimated the time delays, that is, the relative phase between 
the SN images, by fitting the unresolved ground-based flux data with 
the model that includes the flux contributions from the four sets of SN 
lightcurves, F(t, A), each one with its own time of B-band maximum, t 
The fit is constrained by imposing a prior on the image ratios at the 
date of the Keck/NIRC2 observations, reported in Supplementary Table 
3. We allow the lightcurve fit parameters to vary within the ranges 
x =[-3.0, 3.0], c = [-0.3, 0.3], but assume them to be the same for all 
four images. This is an excellent approximation in the absence of appre- 
ciable differential reddening in the lensing galaxy, confirmed by the 
result of the fit, c = -0.01 + 0.01, consistent with no colour excess. 

Galactic extinction was included in the model, adopting 
E(B - V)mw = 0.1558 mag, based on the extinction maps in ref. 40. We 
used the wavelength dependence from the dust from ref. 22 and the 
measured mean value of the total to the Galactic selection extinction 
ratio, Ry= 3.1. From the location of the fitted restframe B-band peak 
luminosity in the SALT2 model for the summed fluxes we inferred the 
time of maximum for SN Zwicky, to = 59808.67 + 0.19, corresponding to 
2022 August 17. The total lensing magnification, y= 24.3 + 2.7, was cal- 
culated after standardizing the SN la fitted B-band peak magnitude with 
the standard SALT2 lightcurve shape and colour correction parameters, 
(a, B) = (0.14, 3.1). Throughout the analysis of SN Zwicky, we adopted a 
flat A cold dark matter cosmological model with Hy = 67.4 km s” Mpc? 
and Qu = 0.315 (ref. 24). 

The uncertainties we report for the time delays account for param- 
eter degeneracies in the fit. Corner plots with posteriors for the relative 
time delays of B, C and D with respect to A are illustrated in Supple- 
mentary Fig. 1. The unresolved data used for this analysis are available 
via WISeREP at https://www.wiserep.org/object/21343. It should be 
emphasized that the lightcurve and time-delay fits were carried out 
completely independently from the lens modelling. 


Spectroscopic follow-up 
The first classification spectrum of SN Zwicky was obtained with 
the integral field unit on the SEDM on 2022 August 21. The data were 
reduced using a custom integral field unit pipeline developed for the 
instrument*’. Flux calibration and correction of telluric bands were 
carried out using a standard star taken at a similar airmass. The details 
of all spectroscopic observations are listed in Supplementary Table 1. 

We obtained three epochs of spectroscopy between 2022 August 
21 and September 11 with ALFOSC on the 2.56 m Nordic Optical Tel- 
escope at the Observatorio del Roque de los Muchachos in La Palma 
(Spain) (see http://www.not.iac.es/instruments/alfosc for further 
information). Observations were taken using grism 4, providing wave- 
length coverage over most of the optical spectral range (typically 
3,700-9,600 A). Reduction and calibration were performed using 
Pypelt version 1.8.14. 

We obtained three epochs of spectroscopy between 2022 August 
22 and October 19 with LRIS on the Keck 110 mtelescope*. All spectra 
were reduced and extracted with LPipe*ć. 


Adaptive optics observations 

Observations with the VLT. HAWK-| observations and data reduction. 
We obtained seven epochs of near infrared in YJHK between 2022 
August 23 and September 30 with HAWK-I*”*’ at the ESO VLT at the 
Paranal Observatory (Chile). All observations, except of the first epoch, 
were performed with the ground-layer adaptive optics offered by the 
GALACSI module® ** to improve the image quality. The first three were 
observed in YJH filters and the next four in YJHK, filters. For the first 
three epochs we exposed for 3 x 60 s in the Y and J bands and 6 x 60 
sin the H band. For the fourth epoch we also observe for 10 x 60 s in 
the K, band. To account for the brightness evolution, we exposed for 
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10 x 100s and 6 x 60s for epochs 5 and 6. For the final epoch we also 
increased H-band exposure times to 10 x 60 s. For the HAWK-I observa- 
tions we used offsets of (—115”, 115”) to place the target on the optimal 
detector chip. 

The data used in our work have been reduced using the HAWK-I 
pipeline version 2.4.11 and the Reflex environment”. The data reduc- 
tion included subtracting bias and flat fielding. The world coordinate 
system was calibrated against stars from Gaia. 


MUSE observations and data reduction. We obtained four epochs 
of integral-field spectroscopy between 2022 August 24 and Septem- 
ber 30 with MUSE™ at the ESO VLT. Each pointing has an approxi- 
mately 1’ x 1’ field of view with spatial sampling of 0.2”per pixel and 
covers the wavelength range from 4,700 to 9,300 A with a spectral 
resolution of 1,800-3,600. All observations were performed with 
the ground-layer adaptive optics offered by the GALACSI module” 
to improve the image quality. The integration time of each epoch 
was 1,800 s. 

The data used in our work have been reduced using the MUSE pipe- 
line version 2.8.7” and the Reflex environment”. The data reduction 
included subtracting bias, flat fielding, wavelength calibration and flux 
calibration against spectrophotometric standard stars. Afterwards, 
we improved the sky subtraction with the Zurich Atmosphere Purge” 
module in the ESO MUSE pipeline. The world coordinate system was 
calibrated against stars from Gaia. 


Laser guide star adaptive optics imaging from Keck. The NIRC2 
J-band observations consist of nine images ina five-point dither pattern 
based ona 2” x 2” grid size, to facilitate sky background subtraction. At 
the first dither location we obtained two images: one co-added 600 s 
exposure and one 200 s exposure. At the second location we obtained 
one 200 s exposure. At each of the third, fourth and fifth dither loca- 
tions, we obtained two 200 s exposures. To correct for flat fielding and 
bias we acquired a set of ten bias frames (flat lamp off) and ten dome flat 
frames (flat lamp on). Sky background and dark current were removed 
as part of the sky subtraction, which utilized a different sky map for 
each dither position, created by median combining the frames from all 
other dither positions, excluding the current dither position. The final 
combined image was created by aligning each dither position to each of 
the others using the centroid of the brightest SN image (image A), and 
median combining. The reduction was carried out using custom Python 
scripts. The NIRC2J-band data (Fig. 2) provides the highest-resolution 
(FWHM 0.086”) image in our dataset of the system where the four SN 
images are visible. 


Modelling of the lens galaxy 

The Keck NIRC2 J-band image was used to model the lens galaxy in 
terms of its 0;, semi-minor to semi-major axis ratio q (or €.=1-q) and 
orientation angle @. The mass profile used in our analysis is a singular 
isothermal ellipsoid”°”: 


O: 


-m (1) 
2V qx +y2/q 


K(X, y) = 


where x corresponds to the convergence (that is the dimensionless 
projected surface mass density) and the coordinates (x, y) are centred 
onthe position of the lens centre and rotated anticlockwise by @. The 
projected mass Minside an isodensity contour of the singular isother- 
mal ellipsoid is given by” 


_ 2 DD p 
M= oe 


(2) 


To calculate D, D,and D,,, we assumed a flat A cold dark matter cosmol- 
ogy with Hy = 67.4 km s” Mpc and Qu = 0.315 (ref. 24). 


In addition to the lens mass model, we included light models 
for the lens galaxy and SN host galaxy in the form of elliptical Sérsic 


profiles: 
1/n 
KR) rop] s [C |l. 6) 


where isthe intensity at the half-light radius R., b, = 1.9992n - 0.3271 


(ref. 57) and 
R=,/x2 +I, (4) 


with q; the axis ratio of the Sérsic profile. The SN images were modelled 
as point sources. We used a Moffat PSF with power index 2.94 and FWHM 
0.091” to model the full image. We simultaneously reconstructed the 
lens mass model, SN images and surface brightness distributions of 
the lens and host galaxy. The lens mass model is constrained only by 
the positions of the lensed SN images and not by their fluxes, since the 
latter may be considerably affected by substructures, such as stars, in 
the lensing galaxy. Our model contains 13 nonlinear free parameters: 
Og, $, q, x,y for the lens mass model, R., n, Ds, Js, Xs, Ys for the lens light 
model and xy, Ysv for the SN position in the source plane. Our results are 
obtained using lenstronomy (https://lenstronomy.readthedocs.io/en/ 
latest/), an open-source Python package that uses forward modelling 
to reconstruct strong gravitational lenses”. The result of the fit and 
comparison with the observations is shown in Supplementary Fig. 2. 
Asacross-check, we also modelled the HST photometry data and redid 
the analysis with LENSMODEL**””, which produced consistent results. 

The resulting best-fit values for the lens mass and light profiles are 
summarized in Supplementary Table 2. Additionally, we derived the 
gravitational mass within the isodensity contour given by the critical 
line of the lens (with radius 6, = 0.167”, corresponding to 0.628 kpc) to 
be M = (7.82 + 0.06) x 10° M,. 

Supplementary Table 3 displays the observed time delays 
and individual fractional flux contributions from each SN image as 
detailed in Lightcurve fit, magnification and time-delay inference 
(təs and fobs), compared with the predictions from the lens model 
(tmoa ANd faoa). Here, time delays are given with respect to image A: for 
example, t;= t;— t,. The observed fractional flux ratios are measured 
from the Keck J-band image after subtracting the lens galaxy light, and 
the model predictions, f,,,4, are computed from the magnifications 
predicted by the lens model, f;= 1/21. In addition to the uncertainties 
obtained from theJ-band image analysis, we make a conservative error 
estimate by also including the scatter in f,,, and f,,,4 obtained when 
modelling the system using data from the HST optical filters F475W, 
F625W and F814W, as well as the Keck near-infrared J-band data. Using 
this approach, we also take into account possible error contributions 
from uncertainties in dust extinction, lens galaxy subtraction and lens 
mass modelling. 

The observed individual image magnifications can be obtained 
from fs by multiplying the individual fractional fluxes by the total 
observed SN magnification of pi" = 24.3 + 2.7. The total lens model 
image magnification, p'°* „is sensitive to the lens mass slope (which is 
unconstrained by the imaging data), such that flatter halos predicts a 
larger magnification. For an isothermal lens, pi" , = 14.9 + 0.9, indicat- 
ing a flatter halo profile (see also ref. 8). However, the predicted flux 
ratios remain approximately constant, which means that to first 
approximation we can multiply the derived fmoa Dy an arbitrary y'" In 
contrast to the predicted individual fractional flux contributions, the 
observed flux is dominated by images A and C. To match the observed 
values, substructure lensing is needed to additionally magnify images 
A and C with factors of >4 and >2, respectively, compared with any 
additional substructure (de)magnifications of images B and D. 

To check the impact of added macro lens model complexity, we 
have studied cases where the lens mass distribution is modelled with 
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two matter components: one where the surface mass density follows 
the lens light distribution (with an arbitrary mass-to-light ratio, pos- 
sibly interpreted as a baryonic mass component), and a second one 
introducing a dark matter halo with additional flexibility on density 
profiles. In spite of the added extra complexity, the quality of the fit 
to the SN image positions does not improve, and induces shifts in the 
predicted flux ratios only below 5%, that is, very small in comparison 
with the observed flux ratio anomalies. 

We do not detect any time dependence in the flux ratios between 
the Keck and HST images, observed just a month after the lightcurve 
peak, 6 d apart, or any other anomalous variation in the unresolved 
~80-d-long lightcurve shown in Fig. 3. Comparing the expected size 
of SN photosphere with the Einstein radii of compact objects in the 
lensing galaxy, we infer that if the discrepancy from the smooth lens 
model is caused by microlensing the deflectors must exceed 0.02 M,. 

With the aim of distinguishing between lensing effects by stars or 
larger substructures, we investigated the upper limit of image splitting 
for the brightest image, A. We approximated the maximum Einstein 
radius of a large structure at the position of image A by putting an 
upper limit on the difference in the FWHM between PSF, and PSFacp. 
Using equation (2), we inferred a 95% confidence upper limit on the 
substructure’s mass within its Einstein radius. Combined with the lower 
mass limit derived from the lack of flux ratio variability, this constrains 
the mass of the substructure deflector in the line of sight to image A to 
0.02 < M/M, <3 x 10”. 

We conclude that, similarly to the case of iPTF16geu, a smooth lens 
density fails to explain the individual image magnitudes and additional 
substructure lensing is needed. Since the observed properties of lens 
systems to first order depend only on the integrated mass within the 
images and/or the surface mass density of the lens at the image posi- 
tions or inthe annulus between the images’, small-image-separation 
systems provide a unique probe of the central regions of gravitational 
lensing galaxies. The probability for lens substructures, such as stars, 
to accommodate the needed (de)magnifications for a range of lens 
density slopes will be investigated in a future lens modelling paper. 


Lens galaxy photometry 

We retrieved science-ready co-added images from Pan-STARRS 
DR1”. Using a circular aperture with a radius of 1.1”, we obtain the 
following apparent magnitudes of the lens galaxy: g = 22.09 + 0.09, 
r=20.71+ 0.02, i= 20.14 + 0.02,z=19.84 + 0.02 and y = 19.63 + 0.05 (all 
errors are of statistical nature). We model the spectral energy distribu- 
tion with the software package Prospector version 1.1. Prospector uses 
the Flexible Stellar Population Synthesis (FSPS) code“ to generate the 
underlying physical model and Python-FSPS™ to interface with FSPS 
in Python. The FSPS code also accounts for the contribution from the 
diffuse gas (for example, H II regions) based on the Cloudy models from 
ref. 65. Furthermore, we assumed a Chabrier initial mass function®® 
and approximated the star-formation history by a linearly increasing 
star-formation history at early times followed by an exponential decline 
at late times (functional form t exp(-t/T)). The model was attenuated 
with the ref. 67 model. 

The best-fit galaxy model points to a moderately massive galaxy 
with stellar mass log M,/M, = 10.1+ 0.3. The other parameters, suchas 
star-formation rate, attenuation and age are poorly constrained and 
we report their values only for the sake of completeness: star-formation 
rate = 157% Mo yr™, E(B—V)4 = 0.6704 mag, age = 1.6*>) Gyr. We 
acknowledge that the aperture might not encircle the entire galaxy 
and, therefore, we might underestimate the stellar mass of the lens. 
Changing the radius of the measurement aperture to 4” increases the 
mass by -0.3 dex. This large aperture also includes the contribution of 
the SN host galaxy and therefore overestimates the stellar mass of 
the lens galaxy. Nonetheless, this upper bound does not alter our 
conclusion about the compact nature of the lensing galaxy, relative to 
other lensing systems. 


Host galaxy photometry 

Numerous studies have shown that the peak absolute magnitudes of 
type la supernovae depend on their host galaxy masses (see, for exam- 
ple, ref. 25). Although the host and the lens galaxy are well separated 
inthe HST images, both galaxies have diffuse emission extending well 
beyond the separation of the two galaxies. To first order, we can remove 
the contribution of the lens by subtracting the lens images predicted 
by lenstronomy from the HST images. Using a circular aperture witha 
1’ radius and appropriate zero points from HST, we measure for the 
host 22.31 + 0.11, 21.81 + 0.06 and 21.13 + 0.06 mag in F475W, F625W 
and F814W, respectively. Fitting the spectral energy distribution with 
Prospector with the same assumptions as in the previous section gives 
a galaxy stellar mass of log M, /Mg = 9.614. A closer inspection of the 
subtracted HST images reveals that the 1” aperture does not encircle 
the total emission of the host. Increasing the aperture radius to 1.5” 
encircles almost the entire host galaxy (F475W = 21.80 + 0.09, 
F652W = 21.38 + 0.04, F814W = 20.74 + 0.03), but it also includes 
non-negligible contribution from residuals of the lens galaxy. The 
galaxy mass increases marginally to log M/M, = 9.7703 but is consistent 
with the previous measurement. This suggests that a mass-step cor- 
rection is not required. To obtain a more robust estimate of the galaxy 
mass, near-infrared observations after the SN has faded are required. 


Data availability 
The reduced spectra and lightcurves used in the paper are avail- 
able through the WISeREP repository® at https://www.wiserep.org/ 
object/21343. The raw VLT data are also available from the ESO Science 
Archive Facility, http://archive.eso.org/eso/eso_archive_main.html, 
program ID:109.234A. 


Code availability 

The data were analysed using the public codes SNTD (https://sntd. 
readthedocs.io/en/latest/) and lenstronomy (https://lenstronomy. 
readthedocs.io/en/latest/). The Keck AO J-band images and analy- 
sis software can be requested from the first author. The Liverpool 
image reduction software will be made available upon publication by 
K-Taggart et al. 
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